Direct extraction of the Eliashberg function for electron-phonon coupling: A case study of Be(lOlO) 
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We propose a systematic procedure to directly extract the Eliasiiberg function for electron-piionon coupling 
from high-resolution angle-resolved photoemission measurement. The procedure is successfully applied to the 
Be(lOlO) surface, providing new insights to electron-phonon coupling at this surface. The method is shown to 
be robust against imperfections in experimental data and suitable for wider applications. 
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Electron-phonon coupling (EPC) is the basis for many in- 
teresting phenomena in condensed matter physics such as con- 
ventional superconductivity. Its possible role in the high- 
Tc cuprates is also being actively discussed [1-3]. Exper- 
imentally, recent advances in high-resolution (both energy 
and momentum) angle-resolved photoemission spectroscopy 
(ARPES) have stimulated many studies on EPC in various 
systems [3-12]. These ARPES measurements usually yield 
the mass enhancement factor A [13], which characterizes the 
strength of EPC, along with some primitive information about 
its spectral structure such as the dominant phonon mode. It 
is highly desirable to take advantage of these high-resolution 
data and obtain full characteristics of EPC. 

Theoretically, the full characteristics of EPC are described 
by the Eliashberg function a^F(a); e, k) [13], the total transi- 
tion probability of a quasi-particle from/to the state (e, k) by 
coupling to boson (phonon) modes of frequency co [14]. Es- 
sentially all physical quantities related to EPC can be deduced 
from the function. For instance, the mass enhancement factor 
A is related to the Eliashberg function by [13], 
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In this Letter, we present a systematic procedure to di- 
rectly extract the Eliashberg function from the high-resolution 
ARPES data. The Maximum Entropy Method (MEM) [15] is 
employed to overcome the numerical instability inherent to 
such efforts [16]. With Eq.l, our procedure also gives a reli- 
able estimate of the mass enhancement factor A, without mak- 
ing ad hoc assumptions on phonon model or requiring low 
temperature measurement [5-12]. 

This procedure is illustrated using new high-resolution 
ARPES data at the Be(lOlO) surface. This system is ideal 
for testing the new procedure because Be has a broad phonon 
band (~ 80meV), thereby removing the need for super-high- 
energy resolution. Also, results of measurements [10, 17] and 
theoretical calculations [18] for this surface have been pub- 
lished, providing references for comparison. 

The photoemission experiments were performed at the Ad- 
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Figure 1: (Color online) (a) Momentum distribution curves (MDCs) 
of SI surface state of Be(lOlO) along AT. The numbers denote the 
initial state energies in meV. The solid lines show the Lorentzian 
fittings, (b) Quasi-particle dispersion determined from the MDCs 
(circles). The solid lines show fittings to the dispersion with the pro- 
cedure detailed in the text using different bare quasi-particle disper- 
sions eo(k). The parameters Vf (eV ■ A/ft) and [3 (eV ■ A^) of eo(k) are 
shown in the inset. The dashed line indicates eo(k) that results in the 
best fit. 



vanced Light Source (ALS) on Beamline 10.0.1 using a dis- 
play high-resolution Scienta 2002 energy analyzer at 24 eV 
photon energy with total energy resolution lOmeV and an- 
gular resolution +0.15° in 6 x 10"" Torr vacuum and at 
T - 30 K. The cleaning procedure for the Be(lOlO) sample 
was described earlier [5]. 

Figure 1(a) shows the momentum distribution curves 
(MDCs) of the Be(lOlO) SI surface state. The quasi-particle 
dispersion e(k) is determined from MDCs by Lorentzian fit- 
tings [8], and is shown in Fig. 1(b). The real part of the self 
energy is calculated by ReS(e, k) = e(k) - eo(k), where eo(k) 
is the bare quasi-particle dispersion, i.e., without EPC. Within 
the small energy scale considered, it is sufficient to approxi- 
mate e()(k) X! -hvpik - kp) + Pik - kp)^. Following the ex- 
trapolation procedure that will be detailed later, we find that 
hvF = 2.95 eV • A and yS = 6.74 eV ■ A^ provide the best fit to 
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Figure 2: (a) Real part of the self energy. Solid line shows the 
MEM fitting. The error bars cr, = 1 meV are determined from 
the noise level of the data, (b) Extracted Eliashberg function. The 
dashed Une shows the constraint function m(w) used: Wd = 15 meV, 
w„ = 100 meV, and mo = 0. 15 (See Eq. 6). The extracted EUashberg 
function is also shown in (a) to indicate the origin of its structures in 
the self-energy (note the x-axis is reversed). 



the data, as shown in Fig. 1(b). The resulting real part of the 
self energy for the SI state is shown in Fig. 2(a). 

The real part of the self energy is related to the Ehashberg 
function by [13]: 

daja^Fiaj;e,k)K(-^,-^j, (2) 

where K(y,y') = dxfix - y)2y'/(x^ - y'^) with fix) be- 
ing the Fermi distribution function. Eq. 2 is valid for a nor- 
mal Fermi Uquid coupled with boson field(s). For a nor- 
mal metal that has no sharp electronic structure around the 
Fermi surface within the small energy scale of the order of 
the phonon bandwidth (80 meV), dependence on the initial 
state energy e can be ignored because da^Ficij; e,k)/d(jj » 
da^F{oj; e,k)/de [13]. Thus, a^F{a); e,k) in Eq. 2 can be re- 
placed by the EUashberg function at the Fermi level a^F(a)), 
defined as a^F(co) = a^Fico; e^, k). 

With this simphfication, Eq. 2 poses an integral inversion 
problem for extracting a^F(a)) from the real part of the self 
energy. The most straightforward way to do the inversion is 
the least-squares method that minimizes the functional: 
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against the Eliashberg function a^F(a>), where D, are the data 
for the real part of the self energy at energy e,-; ReS(e,) is de- 



fined by Eq. 2 and is a functional of a^FicJ); cri are the error 
bars of the data; and No is the total number of data points. 
Unfortunately, such a straightforward approach fails because 
the inversion problem defined by Eq. 2 is unstable mathemat- 
ically and the direct inversion tends to exponentially amplify 
the high-frequency data noise, resulting in unphysical fluctua- 
tions and negative values in the extracted Eliashberg function. 

To overcome the numerical instability in the direct inver- 
sion, one needs to incorporate the physical constraints into 
the fitting process. The most important constraint is that the 
Eliashberg function must be positive. To do this, we em- 
ploy the Maximum Entropy Method (MEM) [15], which min- 
imizes the functional: 
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where is defined in Eq. 3, and S is the generalized 
Shannon- Jaynes entropy. 
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The entropy term imposes physical constraints to the fitting 
and is maximized when a^F(ci)) = mico), where mioS) is the 
constraint function and should reflect our best a-priori knowl- 
edge for the Ehashberg function of the specific system. In this 
study, we use the following generic form: 
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which encodes the basic physical constraints of the Eliashberg 
function: (a) it is positive; (b) it vanishes at the limit w — > 0; 
(c) it vanishes above the maximal phonon frequency. 

The multiplier a in Eq. 4 is a determinative parameter that 
controls how close the fitting should follow the data while not 
violating the physical constraints. When a is small, the fitting 
will follow the data as closely as possible, and when a is large, 
the extracted Ehashberg function will not deviate much from 
m(<jS). There exist a number of schemes (eg. historic, clas- 
sic, Bryan's method, etc.) that choose the optimal value of a 
based on the data and the constraint function m{oS) [15]. In 
this study, the classic method is used. 

To minimize Eq. 4, Eqs. 2 and 5 are discretized using the it- 
erated trapezoid rule for the integral over oj, and a^F{co) is op- 
timized using a Newton algorithm developed by Skilling and 
collaborators that searches along three prescribed directions 
in each iterative step [15]. Detailed reviews of the algorithms 
can be found in Ref. 15. 

Figure 2(b) shows the extracted Eliashberg function of the 
Be(lOlO) SI surface state and the constraint function m{io). 
It resolves a number of peaks approximately at 40, 60, and 
75 meV. Low energy modes with oj < 30meV are also 
evident. When compared with the first-principles calcula- 
tions of the phonon dispersion of this system [18], the re- 
solved peaks correspond well to those surface phonon modes. 
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Figure 3: Various robustness tests to tlie fittings. The red line in 
each figure indicates the nominal result as presented in Fig. 2(b). (a) 
Extracted Eliashberg functions for different values of a. (b) Extracted 
Eliashberg functions for different constraint function parameters. For 
each curve, one of parameters used in Fig. 2(b) is changed: A. ojd => 
0; B. ai,„ => 200 meV; C. niQ => 0.3. The corresponding values of a 
determined by the classic method are: A. a = 1.13; B. a = 0.22; C. 
a = 0.43. The nominal result has a = 1.44. (c) Extracted Eliashberg 
functions of the fittings A-E shown in Fig. 1(b) for different bare 
particle dispersions, (d) Extracted Eliashberg function (black line) 
after dropping the abnormal data point at ~ 47 meV in Fig. 2(a). 



The extracted Eliashberg function automatically cuts off at 
~ 80meV, which is also consistent with the calculation [18]. 
The mass enhancement factor A is calculated by Eq. 1, yield- 
ing a value 0.68 + 0.08, consistent with A = 0.65 obtained 
in previous measurements using temperature-dependent hne 
shapes [10]. 

We have carried out systematic tests to assess the effects of 
various parameters involved in the fitting and the robustness 
of the procedure against data imperfections. These tests are 
summarized in Fig. 3. Based on these tests, which will be 
detailed below, we conclude that the fitting procedure can be 
well controlled to provide reliable physical insights. 

The most important fitting parameter is the multiplier a. 
Figure 3(a) shows the fitting results as a function of a. Chang- 
ing the value of a does not change qualitative features of the 
extracted Eliashberg function such as the number and the po- 
sitions of the peaks, although the quantitative changes of the 
"contrast" are evident. Furthermore, the estimate of A is not 
sensitive to the value of a: for a changing from 0.01 to 100, A 
varies only between 0.64 and 0.69. 

The classic method determines the optimal value of a based 
on the data and the constraint function. Inevitably, the param- 
eters mo, a>R, and w,„ in the constraint function influence the 
decision of classic method, as demonstrated in Fig. 3(b). A 
proper choice of these parameters is important to ensure that 
the algorithm makes the correct decision. In Fig. 2(b), mo 
is roughly the average height of the Eliashberg function, w„, 
is slightly higher than the maximal phonon frequency, and a 
small but nonzero value of a>D is chosen to suppress the arti- 
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Figure 4: (Color online) A test to the data extraction procedure. Solid 
line shows the pre-defined Eliashberg function. The filled circle-line 
shows the extracted Eliashberg function. Dashed line shows m(a>) 
with parameters cjo = 50meV, a)„, = 100 meV, and mo = 0.24. 
Inset shows the data presented to the program and the MEM fitting. 
Gaussian noise with cr, = 1 meV has been added to the self-energy 
data, r = 30 K. 



fact near the zero frequency as seen in in Fig. 3(b) for ud - 0. 
In this way, we have a constraint function that is close enough 
to the real Eliashberg function and is still sufficiently struc- 
tureless for an unbiased fitting. 

The parameters in the bare particle dispersion, Vf and /3, 
are determined from extrapolation to the higher initial state 
energy. Vf and /3 are inter-dependent because given a value of 
/?, there is only one value of Vf (within a small window) that 
can yield the real part of the self energy which has the correct 
asymptotic behavior and can be fitted by Eq. 2. To find the op- 
timal bare quasi-particle dispersion, we tried a number pairs of 
{vf, /3) to generate the real part of the self energy, then ran the 
MEM fitting within e < 150 meV. The optimal pair of {vf, p) 
is the one that provides the best fit to the dispersion data in the 
larger energy window (300 meV). This procedure is demon- 
strated in Fig. 1(b). The coiTesponding extracted Eliashberg 
functions for different values of {vf, /?) are shown in Fig. 3(c). 
It can be seen that fittings with less optimal values of {vf, /?) 
yield result rather close to the nominal one. This removes the 
need for high accuracy in determining the bare quasi-particle 
dispersion. 

The MEM fitting is rather robust against the data imper- 
fections such as accidental data anomaly, as demonstrated in 
Fig. 3(d). The robustness has two origins: (a) physical con- 
straints built in the fitting automatically filter out those un- 
physical data fluctuations; (b) although one data point may be 
anomalous, it is statistically less likely that a number of data 
points become anomalous simultaneously in a coherent way. 

Figure 4 provides a further test on the robustness of the 
MEM fitting. Here, we generate a set of data of the real part of 
the self energy from a pre-defined Eliashberg function. Gaus- 
sian noise is added to simulate the realistic experimental sit- 
uation. Figure 4 shows a typical result of the test. It can be 
seen that the MEM successfully extracts the overall qualita- 
tive features of the predefined Eliashberg function although 
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the data are rather noisy. The mass enhancement factor cal- 
culated from Eq. 1 is 0.35 ± 0.05, which is very close to the 
exact value 0.35. 

Our MEM fitting scheme can be further improved by using 
a better constraint function m(w). For instance, if the phonon 
density of states F(a>) is known from other reliable sources, 
one can use m(w) = af-^ x F(ci>) with af^ being determined 
by optimizing its posterior probability [15]. Such a scheme 
should allow a seamless integration of the existing knowledge 
of the fine spectral structure and the newly extracted coupling 
strengths. 

Our systematic procedure to extract the Eliashberg function 

has a number of advantages, (a) ARPES has much wider ap- 
plicability than the traditional method using the single-particle 
tunneling characteristics [13]. (b) ARPES allows measure- 
ments along different directions, so the Eliashberg function 
a^F{io; e, k) on the whole Fermi surface could be determined, 
(c) Equation 2, the theoretical basis of our method, makes only 
minimal assumption on the nature of the system, i.e., a normal 
Fermi liquid without strong electronic structure near the Fermi 
energy, (d) The procedure only utiUzes the data for ReS(e), 
which is easier to determine and less prone to the imperfec- 
tions than ImX(e). 

Here, our case study also provides new insights to EPC 
at the Be(lOlO) surface, (a) The extracted Eliashberg func- 
tion is significantly different from the simple fictitious mod- 
els (e.g., Debye and Einstein models) previously used to 
interpret the data [8, 10]; (b) More than 75% (0.5 out of 
0.68) of the enhanced A is contributed by the low frequency 
surface modes below 50meV that are not present in the 
bulk phonon spectrum. Similar behavior is also observed in 
Be(OOOl) surface [14, 20]. This raises the question whether 
the low frequency surface phonon modes are responsible for 
the large mass enhancement factors observed in many metal 
surfaces [4—11]; (c) The average phonon frequency [19], 
Inwiog = (2//i) dLij[a-^F(cij)/a)]lna), is calculated to be 
29meV, which is substantially smaller than its bulk value 
(~ 60meV) [5]. This reduces the estimated Tc for possible 
superconductivity. 

In summary, we have proposed a systematic way to ex- 
tract the Eliashberg function from the high-resolution ARPES 
data. The MEM is employed to overcome the data imperfec- 
tions and numerical instability. By using this new technique, 
we have provided new insights to EPC at the Be(lOlO) sur- 
face. We expect the technique to be useful in many situations, 
for instance, in the study of the possible role of EPC in high- 
temperature superconductivity. 
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